clearvars -except Globaloption option

if exist('Globaloption', 'var') == 0
    Globaloption.savefig = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This code reads in the VAR dynamics and produce the Campbell-Shiller decomposition in the paper
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Load data
data = readtable('../../data/main_VAR/US_full.xlsx');
date = data.year;

startdatenum = 1793;
enddatenum = 1946;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
T_post = length(date(startdate:enddate));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option.figure = 0;

Tstart = startdate; 
Tend = enddate; 
T = Tend - Tstart + 1;

date = date(startdate:enddate);

data = data((data.year >= startdatenum - 1), :);
data = data((data.year <= enddatenum), :);

taxrevgdp = data.tau(2:end); % government tax revenue to GDP ratio
taxrevgdp0 = data.tau(1); % government tax revenue to GDP ratio
spendgdp = data.g(2:end); % government spending before interest exdp. to GDP ratio
spendgdp0 = data.g(1); % government spending before interest exp. to GDP ratio
surplusgdp = taxrevgdp - spendgdp; % primary surplus to gdp
rpcgdpgr = data.x(2:end); % real gdp growth (in logs)
inflation = data.pi(2:end); % growth in GDP price deflator (in logs)
ynom1 = data.x1yrCMT(2:end) / 100; % nominal yield on a 1-year Treasury bill, expressed per annum
ynom10 = data.x10yrCMT(2:end) / 100; % nominal yield on a 10-year Treasury bond, expressed per annum

ynom1 = log(ynom1 + 1);
ynom10 = log(ynom10 + 1);

pdm = data.pdm(2:end); % log price-dividend ratio on CRSP vw stock index; dividend is annual and seasonally adjusted
divgrm = data.divgrm(2:end); % nominal log dividend growth; dividend is annual and seasonally adjusted


deltalogg = [log(spendgdp(1)) - log(spendgdp0); log(spendgdp(2:end)) - log(spendgdp(1:end - 1))];
deltalogtau = [log(taxrevgdp(1)) - log(taxrevgdp0); log(taxrevgdp(2:end)) - log(taxrevgdp(1:end - 1))];

yieldspr = ynom10 - ynom1;
y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
yspr0 = mean(yieldspr);
pi0 = mean(inflation);
x0 = mean(rpcgdpgr);
surplus0 = mean(surplusgdp);

% growth of log div to gdp ratio
divgrm = divgrm - inflation;

A0m = nanmean(pdm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[mean(deltalogtau(date <= 1860)) mean(deltalogg(date <= 1860))]
[mean(deltalogtau(date > 1860)) mean(deltalogg(date > 1860))]

tau0 = mean(deltalogtau(date <= 1860));
g0 = mean(deltalogg(date <= 1860));
Tpre = length(date(date <= 1860));

ttspre(1) = taxrevgdp(1); 
gtspre(1) = spendgdp(1); 

for t = 2:Tpre
    gtspre(t) = gtspre(t - 1) * exp(deltalogg(t) - g0);
    ttspre(t) = ttspre(t - 1) * exp(deltalogtau(t) - tau0);
end

tau0 = mean(deltalogtau(date > 1860));
g0 = mean(deltalogg(date > 1860));
Tpost = length(date(date > 1860));

ttspost(1) = taxrevgdp(date == 1861); 
gtspost(1) = spendgdp(date == 1861); 

for t = 2:Tpost
    gtspost(t) = gtspost(t - 1) * exp(deltalogg(Tpre + t) - g0);
    ttspost(t) = ttspost(t - 1) * exp(deltalogtau(Tpre + t) - tau0);
end

gts = [gtspre / exp(mean(log(gtspre))), gtspost / exp(mean(log(gtspost)))];
tts = [ttspre / exp(mean(log(ttspre))), ttspost / exp(mean(log(ttspost)))];

deltalogg = [deltalogg(1:Tpre) - mean(deltalogg(1:Tpre)); deltalogg(Tpre + 1:end) - mean(deltalogg(Tpre + 1:end))];
deltalogtau = [deltalogtau(1:Tpre) - mean(deltalogtau(1:Tpre)); deltalogtau(Tpre + 1:end) - mean(deltalogtau(Tpre + 1:end))];
g0 = 0;
tau0 = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Read in debt/gdp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

gdebt = data.debt_gdp_sargenthall;

plot(date, taxrevgdp, date, spendgdp)
plot(date, tts, date, gts)
[mean(tts(date <= 1860)), mean(tts(date > 1860))]
[mean(gts(date <= 1860)), mean(gts(date > 1860))]


run ../../tools/run_var_estimate;

[mean(log(taxrevgdp(date <= 1860))), mean(log(taxrevgdp(date > 1860)))]
[mean(log(spendgdp(date <= 1860))), mean(log(spendgdp(date > 1860)))]

X2t(:, coint_tax) = [log(taxrevgdp(date <= 1860)) - mean(log(taxrevgdp(date <= 1860))); log(taxrevgdp(date > 1860)) - mean(log(taxrevgdp(date > 1860)))];
X2t(:, coint_spending) = [log(spendgdp(date <= 1860)) - mean(log(spendgdp(date <= 1860))); log(spendgdp(date > 1860)) - mean(log(spendgdp(date > 1860)))];

save(['MAT/res_annual_VAR_nodebt_US_before1946_demean.mat'], '-regexp', '^(?!(option|Globaloption)$).');

ttime = [startdate:enddate];

output_rp = 0.03;


%% Estimation
run ../../tools/cs_estimation.m

%% steady-state upper bound calculation
upper = exp(pxbar) * (mean(taxrevgdp - spendgdp));

gdpreturn = output_rp + y0nom_1 + yspr0 - pi0;


save MAT/USprepara_demean.mat k0x k1x s upper pdX pxbar gdpreturn CFT CFG DR


%% bootstrap
run ../../tools/BS_compute.m;

save MAT/benchmark_US_before1946_demean.mat k0x k1x s upper pdX pxbar gdpreturn s std_coeff gdebt output_rp pdT pdG CFT CFG DR taxrevgdp spendgdp


